diff --git a/sklearn/linear_model/coordinate_descent.py b/sklearn/linear_model/coordinate_descent.py
index 7d65e8038..2f5cb95e2 100644
--- a/sklearn/linear_model/coordinate_descent.py
+++ b/sklearn/linear_model/coordinate_descent.py
@@ -682,7 +682,6 @@ class ElasticNet(LinearModel, RegressorMixin):
 
         Notes
         -----
-
         Coordinate descent is an algorithm that considers each column of
         data at a time hence it will automatically convert the X input
         as a Fortran-contiguous numpy array if necessary.
@@ -690,7 +689,6 @@ class ElasticNet(LinearModel, RegressorMixin):
         To avoid memory re-allocation it is advised to allocate the
         initial data in memory directly using that format.
         """
-
         if self.alpha == 0:
             warnings.warn("With alpha=0, this algorithm does not converge "
                           "well. You are advised to use the LinearRegression "
@@ -709,62 +707,32 @@ class ElasticNet(LinearModel, RegressorMixin):
                              multi_output=True, y_numeric=True)
             y = check_array(y, order='F', copy=False, dtype=X.dtype.type,
                             ensure_2d=False)
+        else:
+            # If check_input is False, ensure X is copied if copy_X is True
+            if self.copy_X:
+                X = X.copy(order='K')
 
         X, y, X_offset, y_offset, X_scale, precompute, Xy = \
             _pre_fit(X, y, None, self.precompute, self.normalize,
-                     self.fit_intercept, copy=False)
+                     self.fit_intercept, copy=True if self.copy_X else False)
         if y.ndim == 1:
             y = y[:, np.newaxis]
-        if Xy is not None and Xy.ndim == 1:
+        if Xy is not None:
             Xy = Xy[:, np.newaxis]
-
         n_samples, n_features = X.shape
         n_targets = y.shape[1]
 
         if self.selection not in ['cyclic', 'random']:
-            raise ValueError("selection should be either random or cyclic.")
+            raise ValueError("selection should be either 'cyclic' or 'random';"
+                             " got (selection=%r)" % self.selection)
 
-        if not self.warm_start or not hasattr(self, "coef_"):
-            coef_ = np.zeros((n_targets, n_features), dtype=X.dtype,
-                             order='F')
-        else:
-            coef_ = self.coef_
-            if coef_.ndim == 1:
-                coef_ = coef_[np.newaxis, :]
-
-        dual_gaps_ = np.zeros(n_targets, dtype=X.dtype)
-        self.n_iter_ = []
+        self.coef_, self.dual_gap_, self.eps_ = map(np.ravel, _path_residuals(
+            X, y, X_offset, y_offset, X_scale, precompute, self.n_alphas,
+            self.alphas, self.l1_ratio, self.eps, self.n_iter, self.tol,
+            self.selection, self.random_state, copy_X=self.copy_X,
+            return_n_iter=True, check_input=False, **params))
 
-        for k in xrange(n_targets):
-            if Xy is not None:
-                this_Xy = Xy[:, k]
-            else:
-                this_Xy = None
-            _, this_coef, this_dual_gap, this_iter = \
-                self.path(X, y[:, k],
-                          l1_ratio=self.l1_ratio, eps=None,
-                          n_alphas=None, alphas=[self.alpha],
-                          precompute=precompute, Xy=this_Xy,
-                          fit_intercept=False, normalize=False, copy_X=True,
-                          verbose=False, tol=self.tol, positive=self.positive,
-                          X_offset=X_offset, X_scale=X_scale, return_n_iter=True,
-                          coef_init=coef_[k], max_iter=self.max_iter,
-                          random_state=self.random_state,
-                          selection=self.selection,
-                          check_input=False)
-            coef_[k] = this_coef[:, 0]
-            dual_gaps_[k] = this_dual_gap[0]
-            self.n_iter_.append(this_iter[0])
-
-        if n_targets == 1:
-            self.n_iter_ = self.n_iter_[0]
-
-        self.coef_, self.dual_gap_ = map(np.squeeze, [coef_, dual_gaps_])
         self._set_intercept(X_offset, y_offset, X_scale)
-
-        # workaround since _set_intercept will cast self.coef_ into X.dtype
-        self.coef_ = np.asarray(self.coef_, dtype=X.dtype)
-
         # return self for chaining fit and predict calls
         return self
 
